library(Seurat)
args=commandArgs(T)

GCMatRData=args[1]
sample=args[2]
workdir=args[3]
regflag=args[4] #T,regress;F,not regress
markerflag=args[5] #T,find marker gene;F,not find
ChoBarFile=args[6] #if not ,fill in "F"
Resol=args[7]  #eg:1

#finish in 2020.10.29

#input:GCMat(RData)
#cluster and find marker gene
#output:
#x.cluster.png,x.BarCluster.xls,x.ClusterTopGene.xls

#Cluster
seuratClust <- function(GCMat,projectName,dims.use,resol,regflag)
{
  
  print("cluster:")
  stri <- paste("dim:",dims.use," resolution:",resol)
  print(stri)
  
  #create
  project.Seurat <- CreateSeuratObject(GCMat,project = projectName)
  
  #normal
  project.Seurat <- Seurat::NormalizeData(object = project.Seurat, normalization.method = "LogNormalize", scale.factor = 10000)

  #MT Gene
  mito.genes <- grep(pattern = "^MT-", x = rownames(x = project.Seurat@assays[["RNA"]]@counts), value = TRUE)
  percent.mito <- Matrix::colSums(project.Seurat@assays[["RNA"]]@counts[mito.genes, ]) / Matrix::colSums(project.Seurat@assays[["RNA"]]@counts)
  project.Seurat <- Seurat::AddMetaData(object=project.Seurat,metadata = percent.mito,col.name = "percent.mito")

  #scale
  if(regflag == "T")
  {
    print("regress:nUMI and percent.mito")
    project.Seurat <- Seurat::ScaleData(object = project.Seurat, vars.to.regress = c("nUMI", "percent.mito"))
  }else
  {
    print("not regress!")
    project.Seurat <- Seurat::ScaleData(object = project.Seurat)
  }
  
  #variable Feature
  project.Seurat <- FindVariableFeatures(object = project.Seurat, 
                                         selection.method = 'mean.var.plot', 
                                         mean.cutoff = c(0.1, Inf), 
                                         dispersion.cutoff = c(0.5, Inf)
  )
  
  #PCA
  #variable genes remove
  VarGenes <- VariableFeatures(object = project.Seurat)
  RemoveGenes1 <- VarGenes[grep("^MT-",VarGenes)]
  RemoveGenes2 <- VarGenes[grep("RPS",VarGenes)]
  RemoveGenes3 <- VarGenes[grep("RPL",VarGenes)]
  RemoveGenesTotal <- c(RemoveGenes1,RemoveGenes2,RemoveGenes3)
  stri <- paste("variable genes:",length(VarGenes)," remove genes:",length(RemoveGenesTotal),sep="")
  print(stri)
  if(length(RemoveGenesTotal) > 0)
  {
    VarGenes1 <- setdiff(VarGenes,RemoveGenesTotal)
    print("remove some variable genes!")
    print(length(VarGenes1))
  }else
  {
    VarGenes1 <- VarGenes
  }
  
  project.Seurat <- RunPCA(object = project.Seurat, 
                           features = VarGenes1, 
                           verbose = FALSE, 
                           npcs = 50)

  #cluster
  project.Seurat <- FindNeighbors(object = project.Seurat, dims = dims.use)
  project.Seurat <- FindClusters(object = project.Seurat, resolution = resol)
  
  #tsne
  project.Seurat <- RunTSNE(object = project.Seurat, dims = dims.use)
  
  return(project.Seurat)
}

#Cluster end.

getMarkerGene <- function(project.Seurat,TopGeneNum,sampleName,workdir)
{
  print("get cluster Marker Gene:")
  stri <- paste("Marker Gene Num:",TopGeneNum,sep="")
  print(stri)
  
  MarkerGeneResult <- FindAllMarkers(project.Seurat)
  TopGeneTableFile <- paste(workdir,"/",sampleName,".ClusterTopGene.RData",sep="")
  save(MarkerGeneResult,file = TopGeneTableFile)
  
  ClusterList <- unique(as.vector(MarkerGeneResult$cluster))
  MarkerGeneTable <- NULL
  for(i in ClusterList)
  {
    ThisClusterTable <- MarkerGeneResult[which(MarkerGeneResult$cluster == i),]
    ThisClusterTopGene <- ThisClusterTable$gene
    GeneNum <- length(ThisClusterTopGene)
    if( GeneNum < TopGeneNum)
    {
      ChoTopGene <- c(ThisClusterTopGene,rep("-",time = TopGeneNum - GeneNum))
    }else
    {
      ChoTopGene <- ThisClusterTopGene[1:TopGeneNum]
    }
    
    MarkerGeneTable <- cbind(MarkerGeneTable,ChoTopGene)
  }
  
  MarkerGeneTable <- as.data.frame(MarkerGeneTable)
  colnames(MarkerGeneTable) <- ClusterList
  TopGeneFile <- paste(workdir,"/",sampleName,".ClusterTopGene.xls",sep="")
  write.table(MarkerGeneTable,TopGeneFile,quote=F,sep="\t",row.names = F,col.names = T)
  
}

#Marker Gene end

#-------------------------------

#main
setwd(workdir)
library(Seurat)
load(GCMatRData)


if(! ChoBarFile == "F")
{
  print("choose barcode")
  ChoBar = readLines(ChoBarFile)
  resolf = as.numeric(Resol)
  GCMat = GCMat[,ChoBar]
}else
{
  print("no choose barcode!")
  resolf = as.numeric(Resol)
}

#cluster
print("cluster:")
all10x.Seurat = seuratClust(GCMat,project=sample,dims.use=1:15,resol = resolf,regflag)
#cluster end.

#save all10x.Seurat
print("save Seurat R data:")
SaveDataFile=paste(sample,".Seurat.RData",sep = "")
print(SaveDataFile)
save(all10x.Seurat,file=SaveDataFile)
#save end.

#plot cluster picture
print("cluster picture:")
CluPicFile=paste(sample,".",Resol,".cluster.png",sep = "")
png(CluPicFile,width=800,height=500)
TSNEPlot(object = all10x.Seurat,label = TRUE,label.size = 6)
dev.off()
#plot cluster picture end.

#Bar - Cluster Table
print("output Barcode-Cluster Table:")
ident <- all10x.Seurat@active.ident
cell.ident <- names(ident)
BarClus <- as.character(ident)
names(BarClus) <- cell.ident

Table <- as.data.frame(BarClus)
Table <- cbind(cell.ident,Table)
colnames(Table) <- c("barcode","cluster")

BarClusFile <- paste(sample,".BarCluster.xls",sep="")
write.table(Table,file=BarClusFile,row.names = F,col.names = T,sep = "\t",quote = F)
#Bar - Cluster Table end.

#marker Gene
if(markerflag == "T")
{
    getMarkerGene(all10x.Seurat,100,sample,workdir)
}else
{
    print("not find marker gene!")
}
#marker Gene end.

